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Final  Report 
on 

Highly  Accurate  Adaptive  Finite  Element  Schemes 
for  Nonlinear  Hyperbolic  Problems 


1  Introduction 


This  document  is  a  final  report  of  research  activities  supported  under  General  Contract 
DAAL03-89-K0120  between  the  Army  Research  Office  and  The  University  of  Texas  at  Austin. 
The  report  describes  work  performed  during  the  period  July  1,  1989  and  June  30,  1992.  The 
Principal  Investigator  of  the  project  was  Professor  J.  T.  Oden.  The  project  supported  several 
Ph.D.  students  over  the  contract  period,  two  of  which  are  scheduled  to  complete  dissertations 
during  the  1992-93  academic  year.  Research  results  produced  during  the  course  of  this  effort 
led  to  six  journal  articles,  five  research  reports,  four  conference  papers  and  presentations, 
one  book  chapter,  and  two  dissertations  (nearing  completion).  More  complete  summaries  of 
these  documents  are  given  later  in  this  report. 


It  is  felt  that  several  significant  advances  were  made  during  the  course  of  this  project  that 
should  have  an  impact  on  the  field  of  numerical  analysis  of  wave  phenomen^  These  include 
the  development  of  high-order,  adaptive,  /infinite  element  methods  for  elasyodvnamic  calcu¬ 
lations  and  high-order  schemes  for  linear  and  nonlinear  hyperbolic  systems.  Also,  a  theory 
of  multi-stage  Taylor- Galerkin  schemes  was  developed  and  implemented7  in  the  analysis  of 
several  wave  propagation  problems,  and  was  configured  within  a  genera/  /jp-adaptive  strat¬ 
egy  for  these  types  of  problems.  Further  details  on  research  results  ar;!d  on  areas  requiring 
additional  study  are  given  in  the  next  section  of  this  report  and  in  an/appendix. 
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2  Review  of  Research  Results 


2.1  General  Goals 

Despite  a  half  century  of  study,  despite  the  introduction  of  modern  computational  techniques 
and  machines,  and  despite  a  multitude  of  papers,  conferences,  and  journals  on  the  subject, 
the  field  of  numerical  analysis  of  complex  wave  phenomena  has  actually  not  progressed  much 
beyond  its  glorious  beginnings  in  the  pre-World  War  II  era  of  computational  modeling.  In 
a  sense,  the  basic  issues  today  are  the  same  as  they  have  been  for  a  half  century:  accuracy, 
stability,  and  consistency  of  the  numerical  approximation  of  the  propogations  of  functions 
which  possibly  possess  discontinuities.  Many  textbooks  are  filled  with  examples  of  successful 
schemes  for  a  simple  one-dimensional  linear  wave  equation,  or  perhaps  Burger’s  equation 
in  one  dimension.  Indeed,  the  theory  of  numerical  analysis  of  wave  phenomena,  as  it  exists 
today,  is  still  basically  a  one-  dimensional  theory,  and  the  most  significant  advances  in  the 
subject  over  the  last  decade  are  the  introduction  of  methods  which  seem  to  work  well  for 
certain  one-dimensional  cases. 

Among  major  goals  of  research  in  this  area  are  the  development  of  highly  accurate,  stable, 
non-oscillatory,  and  convergent  numerical  approximations  to  study  a  multitude  of  features 
of  solutions  of  hyperbolic  systems  of  conservation  laws  that  can  be  used  efficiently  to  model 
wave  phenomena  of  interest  in  science  and  engineering.  But  these  goals  have  proved  to  be 
paradoxical;  e.g.,  monotone  schemes  may  non-oscillatory  and  stable,  but  they  are  only  first- 
order  accurate,  and  higher  order  schemes,  while  providing  higher  accuracy,  are  almost  always 
oscillatory  and  frequently  unstable.  Moreover,  for  nonlinear  hyperbolic  systems,  questions 
of  uniqueness  of  solutions  arise  and  the  additional  requirement  that  the  numerical  solution 
be  “physically  meaningful”  must  be  added  to  the  list  of  criteria. 

The  emergence  of  adaptive  computational  schemes  over  the  last  decade  has  provided  a 
possible  basis  for  achieving  the  historical  goals  in  the  numerical  analysis  of  wave  phenomena: 
manipulate  mesh  parameters  in  such  a  way  that  both  high-order  accuracy  and  stability  can 
somehow  be  achieved.  In  this  way,  the  computational  process  is  optimized,  and  any  tendency 
of  the  scheme  to  oscillate  or  to  lose  accuracy  would,  in  theory,  be  compensated  by  appropriate 
adaptation  of  the  control  parameters.  A  key  to  this  approach  is  the  basic  idea  of  developing 
reliable  a  posteriori  estimates  of  the  numerical  error  in  a  finite  element  approximation,  mesh 
parameters  such  as  the  mesh  size  h ,  spectral  order  of  approximation  p ,  or  the  location  or 
the  relocation  of  nodes,  can  be  adapted  to  keep  error  within  preset  tolerances.  Importantly, 
the  computational  cost  can  also  be  added  to  the  cost  functional  so  that,  at  least  in  theory, 
highly  efficient  schemes  can  be  used  to  achieve  these  classical  goals. 

The  project  summarized  in  this  final  report  had  as  its  basic  objective  the  exploration  of 
new  types  of  /)p-finite  element  methods  for  the  analysis  of  wave  propagation  problems,  with 
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particular  emphasis  on  stress  waves  in  solids.  The  starting  points  for  all  of  the  studies  in 
this  effort  were  the  following:  a  high-level  ftp  data  structure  was  constructed  in  which  the 
mesh  size  ft  and  the  spectral  order  p  within  a  finite  element  mesh  were  designed  to  be  free 
parameters.  The  capability  for  either  using  ft-refinement,  p-enrichment,  or  ftp-refinements 
was  embodied  in  the  data  structure,  the  data  structure  was  designed  to  admit  fairly  general 
a  posteriori  error  estimates,  so  that  some  research  and  experimentation  could  be  done  on 
the  calculation  of  errors  to  drive  any  adaptive  process.  Finally,  the  specifics  of  the  adaptive 
process  were  left  open  so  that  a  number  of  different  strategies  could  be  explored.  With 
regard  to  the  general  approaches  considered  in  the  research,  two  were  explored: 

1.  Use  of  discontinuous  ftp-finite  element  methods  for  nonlinear  hyperbolic  systems,  and 

2.  the  use  of  continuous  ftp- approximations  with  very  high-order  implicit  schemes  for  the 
particular  classes  of  hyperbolic  systems  that  arise  in  the  study  of  stress  waves  in  solids. 

These  two  approaches  are  discussed  in  more  detail  in  the  following  sections  and  some  tech¬ 
nical  details  are  given  in  the  appendices. 

2.2  Stress  Waves  in  Solids:  RK  and  TG  Schemes 

We  briefly  outline  here  the  approach  and  some  results  of  the  research  on  high-order  schemes 
for  calculation  of  transient  phenomena  in  elastic  solids. 

To  construct  a  successful  adaptive  scheme  for  the  hyperbolic  partial  differential  equa¬ 
tions  that  occur  in  linear  elastodynamics,  the  basic  ingredients  mentioned  above  must  first 
be  addressed:  1),  a  functional  ftp-adaptive  data  structures  in  hand,  a  high-order  temporal 
scheme  must  be  identified  to  advance  the  solution  in  time;  2),  an  efficient  method  of  a  pos¬ 
teriori  error  estimation  must  be  developed  to  provide  the  data  for  controlling  the  numerical 
process;  3),  an  adaptive  strategy  must  be  developed  to  control  (optimize)  the  mesh  during 
the  evolution  of  the  wave  phenomena. 

In  the  early  months  of  this  phase  of  the  project,  considerable  effort  was  spent  exploring 
classical  high-order  methods  for  time  integration.  Surprisingly,  relatively  few  existing  meth¬ 
ods  survive  our  criteria  for  efficiency  and  applicability  to  very  large  systems.  The  classical 
Adams-Bashforth  methods,  for  example,  require  a  significant  amount  of  memory  and  are 
only  conditionally  stable.  The  Runge-Kutta  methods,  however,  do  seem  to  offer  a  number 
of  advantages.  They  are  implemented  locally  over  a  time  step,  which  was  convenient  for  ftp- 
adaptive  schemes,  and  they  could  be  used  to  produce  results  of  arbitrary  order  in  time.  A 
great  deal  of  time  was  spent  during  early  phases  of  this  project  encoding  and  experimenting 
various  forms  of  Runge-Kutta  methods,  including  the  singly-implicit  Runge-Kutta  schemes, 
implicit  Runge-Kutta  schemes,  and  traditional  semi-implicit  Runge-Kutta  schemes.  These 
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proved  to  be  somewhat  effective,  but  in  many  cases  also  were  accompanied  by  unpleasant 
oscillations  for  very  high  order  schemes. 

A  new  study  was  initiated  on  a  completely  new  family  of  methods  which  are  of  a  form 
similar  to  so-called  Taylor-Galerkin  methods  used  in  certain  flow  calculations.  Traditionally, 
Taylor-Galerkin  methods  no  higher  than  third-order  appear  in  the  literature,  and  these  are 
known  to  be  very  inefficient  and  expensive.  However,  it  was  observed  that  by  reducing  the 
elastodynamics  problem  to  a  first-order  hyperbolic  system  and  then  applying  the  Taylor- 
Galerkin  strategy,  a  recurrence  formula  could  be  derived  which  could  produce  very  robust 
and  stable  schemes  of  arbitrary  high  order.  Thus  was  invented  the  first  high-order  Taylor- 
Galerkin  schemes  for  wave  propagation.  These  schemes  are  still  under  study,  but  they  are 
known  to  possess  many  attractive  features  and  have  proved  to  be  quite  superior  to  traditional 
Runge-Kutta  methods  in  a  number  of  numerical  experiments. 

In  the  final  months  of  this  phase  of  the  project,  still  another  version  of  the  implicit,  high- 
order,  multi-staged  Taylor-Galerkin  scheme  were  developed  in  which  second-order  hyperbolic 
systems  are  derived  which  are  equivalent  to  the  equations  of  linear  elastodynamics.  It  v. 
observed  that  the  governing  operators  naturally  split  into  similar  component  parts  which 
made  application  of  the  TG  ideas  straightforward. 

Considerable  time  was  spent  on  a  posteriori  error  estimation  for  these  types  of  schemes. 
These  a  posteriori  estimates  were  developed  for  not  only  the  time-dependent  case  but  also 
for  the  elliptic  systems  obtained  during  each  step  of  the  Taylor-Galerkin  approximation. 
The  situation  is  this:  rigorous  mathematical  theory  for  a  posteriori  error  estimates  for  linear 
elliptic  systems  was  developed  and  applied  step-wise  to  the  Taylor-Galerkin  scheme.  In  the 
Taylor-Galerkin  scheme,  the  elliptic  step,  of  course,  involves  the  time  step  so  that  the  a 
posteriori  estimate  does  include  a  description  of  this  mesh  parameter.  A  number  of  oppor¬ 
tunities  for  measuring  temporal  error  also  present  themselves.  In  particular,  the  usual  use  of 
predictor  and  a  corrector  in  time  allow  for  a  fairly  straightforward  estimate  of  the  temporal 
component  of  the  approximation  error.  To  date,  a  number  of  test  problems  have  been  run 
and  results  suggest  that  the  total  error  in  the  two-dimensional  elastodynamics  problem  can 
be  kept  under  control  and  that  fairly  accurate  error  estimates  can  be  obtained. 

Filially,  it  is  necessary  to  address  the  /ip-adaptive  strategy.  This  is  an  area  in  which  con¬ 
siderable  additional  work  reamins  to  be  done.  During  the  course  of  the  project  reported  here, 
a  number  of  so-called  three-step  schemes  were  explored  in  which  criteria  were  established 
for  producing  first  an  /(-refinement  of  a  mesh  and  then  a  p-ennchment  to  control  the  error. 
Various  versions  of  this  strategy  have  been  studied  experimentally.  The  meshes  obtained  by 
these  schemes  are  certainly  not  optimal,  and  some  are  quite  far  from  optimal,  but  they  can 
be  implemented  with’ great  speed  and  thus  the  overall  computational  time  of  implementation 
often  proves  to  be  quite  acceptable. 
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The  status  of  these  schemes  is  this:  a  number  of  papers  on  the  mathematical  ideas  and 
the  implementation  have  been  published,  a  working  research  code  for  two-dimensional  cases 
has  been  developed,  and  a  number  of  test  problems  have  been  run.  Further  work  remains  to 
be  done  on  error  estimation,  adaptive  strategy,  and  on  numerous  details  of  implementation, 
such  as  the  possibility  of  using  domain  decomposition  and  parallelization  during  the  compu¬ 
tational  process.  Mathematical  issues  deserving  study  include  the  study  of  stability  of  the 
TG  schemes,  a  priori  error  estimation,  proof  of  convergence  of  the  /ip-adaptive  strategics, 
and  further  work  on  rigorous  a  posteriori  estimates. 

2.3  Discontinuous  hp  Methods  for  Nonlinear  Conservation  Laws 

As  noted  earlier,  much  as  been  said  in  recent  literature  on  the  numerical  solution  of  hy¬ 
perbolic  conservation  laws  with  regard  to  the  use  of  flux  limiting  methods  and  high-order 
approximations.  This  has  led  to  the  notions  of  TVD  schemes,  ENO  schemes,  etc.,  all  of  which 
seem  to  work  very  well  for  one-dimensional  problems  provided  that  no  boundary  conditions 
of  significance  are  imposed.  To  extend  these  ideas  to  reasonable  two-  and  three-dimensional 
cases  involves  some  significant  generalizations  of  the  existing  theory. 

One  generalization  that  was  felt  deserved  some  study  was  the  notion  of  TVB  (total  vari¬ 
ation  bounded)  algorithms  developed  Cockburn  and  Shu.  While  their  work  itself  originally 
focused  on  one-dimensional  cases,  much  of  their  theory  is  general  and  potentially  extendable 
to  higher  dimensional  cases. 

As  another  thrust  in  the  present  research  effort,  the  study  of  high-order  hp  schemes 
for  discontinuous  Galerkin  approximations  of  nonlinear  hyperbolic  systems  was  undertaken. 
These  were  based  on  generalizations  of  the  Cockburn-Shu  TVB  schemes,  and  employed  high- 
order  hp  quadrilateral  finite  element  approximations  in  two  dimensions.  The  key  to  these 
types  of  methods  is  the  construction  of  a  special  projection  operator  which  maintains  the 
TVB  character  of  the  numerical  solution.  During  this  project,  an  algorithm  for  producing 
such  a  projection  was  indeed  developed  and  has  been  applied  successfully  to  discontinuous  hp 
approximations  of  hyperbolic  conservation  laws.  The  status  of  this  work  is  that  one  paper 
on  the  subject  has  been  published,  another  is  in  preparation.  Several  important  details 
remain  unresolved.  These  include  the  development  of  a  useful  adaptive  strategy  for  these 
techniques,  and  a  numerical  strategy  that  truly  exploits  the  spectral  character  of  these  types 
of  approximations  while  using  the  A-adaptivity  to  control  oscillations  near  shocks  and  contact 
discontinuities. 
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APPENDIX 


High-Order  Taylor-Galerkin  and  Adaptive 
hp-  Methods  for  2nd  order  Hyperbolic 

Systems 


In  this  Appendix  a  brief  summary  of  some  of  the  details  of  approach  and  numerical  results 
of  the  hp-adaptive  Taylor-Galerkin  algorithms  described  in  the  text  is  given.  This  appendix 
is  excepted  from  an  article  on  the  subject  which  is  entitled  "High-Order  Taylor-Galerkin 
and  Adaptive  h-p  Methods  for  Hyperbolic  Systems"  by  A.  Safjan  and  J.  T.  Oden. 


1  INTRODUCTION 


’’  he  use  of  high-order  adaptive  finite  element  methods  for  elliptic  boundary-value 
problems  in  which  the  mesh  size  h  and  the  local  spectral  order/?  of  the  approximation  are 
varied  in  order  to  control  the  approximation  error  have  been  shown  to  produce  exponential 
rates  of  convergence  (e.g.,  [7,  8]).  These  methods  treat  the  mesh  variables  h  and  p  as 
arbitral'  parameters,  and,  through  a  fairly  elaborate  data  structure,  distribute  mesh  sizes 
and  orders  nonuniformly  over  a  mesh  to  control  local  errors,  which  are  estimated  using  a- 
posteriori  estimation  techniques  [7].  Thus,  they  represent  optimal-control  strategies  which 
attempt  to  configure  the  mesh  to  optimize  the  computational  process.  Such  strategies  have 


proved  to  be  effective  for  several  classes  of  problems,  including  problems  in  gas  dynamics 
modeled  by  the  Euler  equations  [1]  and  pioblems  in  compressible  and  incompressible  flow 
(e.g.,  [71). 

For  time-dependent  problems,  however,  the  use  of  high-order  spatial  approximations 
requires  the  use  of  a  balanced  temporal  approximations  that  is  high  order  as  well. The 
search  for  robust  high-order  schemes  that  function  efficiently  on  spatially-nonuniform 
meshes  has  been  an  elusive  one,  and  few  of  the  traditional  schemes  of  high  order  (e.g., 
implicit  Runge-Kutta  methods)  have  proven  to  be  effective  for  these  types  of 
approximations  [9]. 

The  so-called  Taylor-Galerkin  schemes  represent  generalizations  of  the  Lax-Wendroff 
algorithm  and  have  been  used  effectively  for  producing  second-  and  third-order  temporal 
approximations  [5,  3].  However,  no  procedures  forextending  these  techniques  to  temporal 
approximations  of  arbitrary  order  appear  to  be  available. 

In  the  present  work,  we  present  a  new  family  of  stable  high  order  Taylor  Galerkin 
(TG)  methods  for  the  numerical  solutions  of  second-order  hyperbolic  systems.  It  is  shown 
that  for  second  order  systems,  a  mulu-stage  process  can  be  used  that  produces  schemes  of 
order  2s  for  ^-stages,  each  of  which  involves  the  solution  of  a  second-order  system  of 
elliptic  equations.  A  detailed  stability  analysis  is  provided  for  the  case  of  linear  systems 
which  establishes  choices  of  parameters  that  result  in  unconditionally  stable  schemes. 

An  error  estimation  procedure  is  also  presented  which  leads  to  estimates  of  both  the 
spatial  and  temporal  approximation  error.  In  addition,  an  adaptive  algorithm  is  developed 
which  employs  an  ftp- adaptive  Finite  element  method  for  controlling  the  spatial  errors. 

For  focus,  applications  to  linear  elastodynamics  problems  in  two  space  dimensions 
are  described.  The  results  of  numerical  experimnts  on  representative  two-dimensional 
stress  wave  propagation  problet.  are  also  given.  The  results  indicate  that  the  algorithms 
are  capable  of  delivering  high  accuracy  on  meshes  with  very  high  order  spatial 
approximations.  Very  little  oscilations  of  solutions  in  the  vicinity  of  wave  fronts  is 
observed,  despite  the  very  high-order  approximations  and  minimal  numerical  dissipation. 
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2.  Model  2nd  order  Hyperbolic  System:  Equations 
of  Linear  Elastodynamics 


As  a  stoning  point ,  we  consider  equations  of  linear  elasticity  in  the  following  form: 
divT(u)  +  p0b  =  p0dju 

T(u)  =  2pe(u)  +  A7>(e)/  in  £2  (2.1) 

2e(u)  -  Vu  +  VuT 


where : 


Q  is  a 
u 

T(u) 
£(«) 
b(x,  t) 

Po 
p,  A 


domain  in  RN,  N  =  2,3 

=  u(x,  t )  is  the  displacement  vector  at  particle  x  e  12  at  time  l  >  0 

=  the  stress  tensor 

=  the  strain  tensor 

=  the  body  force  field 

=  mass  density 

=  Lame  coefficients 


We  reformulate  equations  (2.1)  in  terms  of  u  only  and  arrive  at 

2 

Podtu  -  DtCDu  =  pnb 
where  D  is  a  generalized  gradient  operator  : 


(2.2) 
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0  0 


and  C  is  a  6  x  6  symmetric  positive  definite  matrix  of  elastic  constants  : 

^1  111  ^1122  ^  1 133  ^1123  ^1131  ^1112 
C2222  ^2233  ^2223  ^2231  ^2212 

C3333  C3323  c333]  c33i2 

^2323  ^2331  ^2312 
^3131  C3U2 

*  «  •  •  t 

^  1212 

For  isotropic  material  elasticities  C-,ju  are  given  in  term  of  Lam£  constants  X  and  ji  by 

Ci}id  =  Sij  8/a  X  +  /a(  8ik  8ji  +  8;i  8jk )  (2.4) 

Equations  (2.2)  are  to  be  solved  in  a  domain  Qcj?w,  N  =  2,  3.  Typically,  two  particular 
cases  are  of  interest: 

•  interior  problems  when  Q  is  bounded 

•  exterior  problems  when  £2  is  a  complement  of  a  bounded  set 

The  initial  boundary  value  problem  is  further  specified  by  introducing  boundary 
conditions.  We  consider  the  following  kinds  of  boundary  conditions: 

1 .  Kinematic  boundary  condition 

u  =  u  on  Tu  (2.5) 

where  u  is  a  prescribed  displacement  vector  on  the  Tu  -  part  of  the  boundary  . 

2.  Traction  boundary  condition 
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T(u)  •  n  =  t 


on  r, 


where  n  is  the  unit  outward  normal  to  the  boundary  and  /  is  a  prescribed  traction  on  the  Tf 
-  part  of  the  boundary  (<90  =  Tu  vj  r>  T,  =  0). 

The  initial  boundary  value  problem  is  completed  by  specifying  initial  conditions  of  the 
form 

u  =  Uq  and  dtu  =  v0  at  t  =  0  (2.7) 


In  the  case  of  homogeneous  boundary  conditions,  the  problem  can  be  cast  into  the 
Hilbert  space  formulation  as  follows  (see  [4]  for  a  detailed  discussion). 

We  introduce 


•  The  Hilbert  space 


//  = 


{u,v)H  =  (u,Mv)(L  2yv 


with  the  weighting  matrix  M  =  (mu  ) 


mkl  ~  Po^ki 


•  Operator  Iff  :  H  z>  D{ Iff) — *  H 


Su  =  -M'DtCDu 


(2.10) 


where  D{8)  is  the  domain  of  8  defined  as 


D(8)  =  {  #/J(Q)  |  DtCDu  e  H } 


(2.14) 


for  the  Dirich  let  boundary  value  problem,  and 


D{8)  =  (  ll\Q)  !  CDu  e  D*(Q)} 


(2.15) 


for  the  Neumann  boundary  value  problem,  where 

/)*(n)  =  {  ue  (L^))M  |  e  H(Q),  Vv  e  n\Q),  (u,  Dv)^m  =  (>•'.  v)(L2yv  } 

(2.16) 

and  M  =  N(N+ 1)/2  .  Note  that  the  boundary  condition  on  u  is  satisfied  in  the  sense  of  the 
trace  theorem,  whereas  the  boundary  condition  on  /  is  interpreted  in  the  sense  of  the 
generalized  Green’s  formula.  For  these  reasons,  the  displacement  boundary  condition  is 
classified  as  the  Dirichlet  boundary  condition  and  the  traction  boundary  condition  as  the 
Neumann  boundary  condition  for  operator  8. 

Within  the  Hilbert  space  formalism,  the  initial  boundary  value  problem  of  linear 
elasticity  can  be  reinterpreted  as  an  abstract  2nd  order  Cauchy  problem  : 


d2u 
dt 2 


8u 


0 


t  >  0 


u 


du 

uo  ♦  df  =  vo 


i  =  0 


(2.17) 


An  //  -  valued  function  of  time  u  =  u(t) 

[0,  oo)3/ — »  u(t)e  ll  (2.18) 

is  called  a  weak  solution  of  (2.17)  if : 

(i)  Uq  and  satisfy  regularity  assumption 

«o^0  €  ff  (2.19) 

(ii)  u(r)  satisfies  regularity  assumption 

ue  C([0,  oo);//)  (2.20) 

(iii)  u(t)  satisfies  (2.21) 

Jo  "r(0./r+  %0  +  (^^0-’))(L2f 
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0 


(2.21) 


-  (v0,  <**(), •  ))(j2yv  = 


for  every  test  function 


0  e  C0(i?  ;  D(ff)  )  n  C2(i?  ;  //  ) 


(2.22) 


Notice  that  this  definition  admits,  in  particular,  solutions  in  the  d'Alembert  sense. 


We  record  now  some  fundamental  results  concerning  operator  iff  and  the  existence  and 
uniqueness  of  weak  solutions  u.  We  restrict  ourselves  to  the  case  of  a  bounded  domain  Q. 

1 .  Operator  8 is  self-adjoint. 


2.  Operator  #1/2  exists. 

3.  The  spectrum  of  til,  of  til),  consists  of  eigenvalues  only.  For  the  Dirichlet  problem  all 
the  eigenvalues  are  positive 


o(2f)  =  {t»,  ,  (02,  ...  } 


0  <  ct>j  <(o2<  ...  <  (on 


and  for  the  Neumann  problem  the  eigenvalues  are  non-negative 


(2.22) 


of  Iff)  =  {  o>0,  (Oj  ,  co2  ,  ...  } 

0  =  (Oq  <  Ct>,  <  col  <  ...  <  (On  - >  oo 


(2.23) 


All  eigenvalues  are  of  finite  multiplicity  and  corresponding  eigenspaces  («n)  are 
orthogonal. 

4.  For  the  Neumann  problem,  the  eigenspace  corresponding  to  zero  eigenvalue,  i.e.,  the 
null  space  of  operator  til,  is  spanned  by  constant  vectors  and  by 


(2.24) 
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5.  Operator  S admits  a  clasical  spectral  decomposition 


(2.25) 

(2.26) 


where  dP n  is  the  orthogonal  projection  on  the  eigenspace  {«rt}  corresponding  to  (on 

dPn{')  =  (*,  «>„  (2.27) 

and  «o  =  1, 0  for  the  Dirichlet  and  Neumann  problems,  respectively. 

5.  A  weak  solution  u  exists  and  is  unique.  Moreover,  it  is  of  the  form 

u(t )  =  cos  ( #1/2 t)u0  -  sin  (lff1/2  t)v0  (2.28) 


6.  If  the  initial  condition  functions  w0  and  vQ  satisfy  an  additional  regularity  assumption 

u0E  D(8M).  v0G  //  (2.29) 

and  the  solution  uE  C1  ( [0,  oo)  ;  H  )  n  C(  [0,  oo) ;  D{  ),  then  the  weak  solution  is 
also  a  solution  with  finite  energy: 

£(•)  =  ll«  ,  \\2h  +  llffi/2 u  W],  (2.30) 
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3  High-Order  Taylor-Galerkin  Methods 


Given  a  bounded  domain  ft  in  RN ,  N  =  2,  3,  we  consider  a  system  of  conservation 
laws  of  the  form 

u tt  +  Fk(Vu)  k  =  0  xe  £2  r>0  (3.1) 

k  =  1,  ....  N 

where  u  =  m(x,  t )  is  a  column  vector  of  M  unknowns,  Fk ,  k  =  I, N  are  vector-valued 
functions  of  Vu ,  commas  denote  the  differentiation  with  respect  to  time  t  and  spatial 
variables  xk  and  the  usual  summation  convention  holds. 

This  system  of  equations  is  accompanied  by  an  initial  condition, 

u(x,  0)=  Uq(jt),  u  t(x,  0)=  v0(x)  xe  Q  (3.2) 

and  by  appropriate  boundary  conditions. 

By  introducing  velocity  v  ~  u  (  as  an  auxiliary  variable,  equations  (3.1)  can  be 
converted  into  the  following  first  order  system  (in  time): 


{v  t  +  F  k(  Vu )  k  ~  0 
-  "  =  0 


xe  Q  r>0  (3.3) 


Finally,  (3.3)  and  (3.2)  is  cast  in  the  form  of  an  abstract  Cauchy  Problem 


r 

U  ,  +  AU  =  0  t  >  0 

U  =  u0  t =  0 


(3.4) 


where 


Fk(V(’))k 

0 


(3.5) 
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'V 

V  is  a  group  variable,  U  = 

|  ,  and  U0  specifies  initial  conditions,  U0  = 

W 

un 
\  °J 

3.1  Taylor-Galerkin  Schemes  for  Nonlinear  Systems 

Given  the  solution  Un  =  U(tn )  at  time  tn  =  nAt,  we  seek  the  next  time  step  solution 
(//)+!  _  U(tn  +  At)  in  the  following  form 


Z  t  - 

XAt2  Z  \  tU  =  Un 

+  M\oAt  unt 

+ 

v'loAr  2U]t 

z2  - 

aa  t2Z2M  =  un 

+  MioAt  Vnt 

+ 

V20A  t2Un{t 

• 

+  M2\At  Zu 

+ 

V2jAr  2Zut 

Zs  - 

AA  t2ZStU  = 

+  MsoAr  Un{ 

+ 

v'soA 1 2Vnn 

+  Ms\AiZ\'t 

+ 

vslA  t2Z\Jt 

+  ... 

+  Ms  1  At  Zs. 1  ,r +■  Vsj- 1  At  “Zs. 1  ,tt 


(3.6) 


where 


Zt  ~  U(tn  +  ci  At),  /  =  1,2, ...,  s  -  1 ,  are  intermediate  solutions  called  “internal 
approximations” 

2.v  =  U(tn  +  At)  =  Un+]  is  the  next  time  step  solution  (i.e.,  cs  =  1) 
t*ij  '  Vij  ’  ^.0  ’  ^(0  ’  ci  e  (=1,2,  ...,  5  ,j=  1,  2, ...,  i  —  1 
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A  e  R+  is  a  stability  parameter 
s  =  number  of  stages 

Coefficients  fJ.-  ,  v -- ,  ni0 ,  v(0  ,  ciy  are  to  be  chosen  so  as  to  obtain  the  highest  possible 
order  of  accuracy,  subject  to  stability  or  other  constraints.  A  free  parameter  A  is  to  be 
chosen  from  stability  considerations. 

It  is  convenient  to  rewrite  (3.6)  in  the  following  compact  form: 


Matrices  v,  fi  and  k together  with  vector  c  completely  characterize  difference  scheme  (3.6). 

A  distinct  feature  of  (3.6)  is  that  coefficient  matrices  v'and^  are  lower  triangular 
matrices  which  makes  the  resulting  scheme  semi-implicit  (i.e,  to  compute  Z,-  it  is  necessary 
to  know  Z/.i,  Zi-2,  ....  Zj,  but  it  is  not  necessary  to  now  ZlJr\).  Moreover,  all  diagonal 
elements  of  v  are  equal  (u,-,-  =  A,  i=l,  2, ...,  s )  and  so  are  those  of  ji  (jin  =  0,  /=  1,  2,  ..., 
5).  This  makes  the  operator  defining  the  left-hand  side  of  each  stage  of  (3.6)  identical  for 
linear  (or  linearized)  problems,  and,  hence,  significantly  reduces  the  cost  of  the  method.  A 
particular  choice  of  zero  diagonal  elements  of  fi  is  made  with  an  eye  on  a  well-posedness  of 
a  typical  one  stage  problem  and  a  possible  splitting  of  the  operator  defining  the  left  hand 
side  of  each  stage. 

To  make  the  t'-th  stage  solution  Z,  m-th  order  accurate,  it  is  necessary  to  satisfy  the 
order  conditions  for  Z;  and  to  make  the  previus  stage  solutions  Z,.|,  Z,-. 2, ....  Zj,  to  be  at 
least  of  the  order  m-\.  (Otherwise,  some  coefficients  have  to  be  set  to  zero,  e.g.,  if  Z,\\  is 
of  the  order  m- 2,  then,  necessarily  Wj.\  -  0).  The  order  conditions  for  Z,  are  obtained  by 
expanding  it  in  Taylor  series  at  Un  : 


Zi  =  UH+(CiAt)U*  +^(ctAt) 2(J*n  +  ...  +  ~  (ciAt)m~f^rUn  +0(Atm+t)  (3.9) 


plugging  (3.9)  into  the  left-hand  side  of  the  z-th  stage  equation 
z,  -  £>=1  (WjAf  ZU  +  vi;At2Z;- (1  )  =  Un  +  Ka&tU*t  +  KiT,Ai2Unu  (3.10) 


and  equating  coefficients  of  powers  of  At  to  zero.  This  leads  to  the  following  system  of 
nonlinear  algebraic  equations 


c?  -  c‘  W  -  k(k-'f£UcJ  lv‘i 


k  ~  1,2,  ...,  m 


Mio  ,  i  =  l 

‘  2  u/o  ,  /  =  2 

.  0.  otherwise 


(3.11) 


Equations  (3. 1 1)  are  referred  to  herein  as  the  order  conditions. 

We  now  proceed  to  derive  coefficients  for  some  particular  schemes.  In  the  sequel  we 
adopt  the  following  notation:  TG(.v,  m)  =  .r-stage  m-th  order  scheme. 
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We  make  Z\  and  Z2  to  be  0(Ap)  and  0(A/4),  respectively,  which  leads  to  the  following 
system  of  equations: 


/iio  =  ci 
2  v'i  0  =  -  2  2, 

62.=  c\ 

^20+^21  =  1 

2(v20+  V21)  +  2^21^1=  1-  21 
6 V2 1  c  1  +  3/^21  =  I-  62. 

12v21cf  +  4ji2i  c\  =  1  -  122. 


The  solution  of  (3.12)  reads: 


C 1  = 

Hio 

=  (6X)1*2 

V10  = 

22. 

1  -  62.  1  1  -  122. 

^21  ~ 

62. 

2  (6A)3/2 

1  1 

-  122.  1  1  -  122. 

V21  = 

4 

62.  ’  3  (62.) 1/2 

1  -  62.  1  1  -  122. 

^20  = 

1  -  - 

61  '  2  (6  2.) 3/2 

1 

,  1  1-122. 

V20  = 

2  ' 

2.  +  -7  -  - 

4  6A 

2  1  -  62. 
3  {61)1'1 


(3.12) 


(3.13) 


A  different  4th  order  scheme  can  be  obtained  by  discarding  the  3rd  equation  of  (3.12) 
(i.e.,  making  Z\  to  be  only  0(A/2)  )  and  by  taking  /i2i  =  0-  The  coefficients  for  this 
scheme  are  given  below 


c\ 


V10  = 
^21  = 


Vio  = 


1  1  -  122. 

2  1  -  62. 

1  -  122.\2 


1  -  62. 
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V2i  = 

1 

3 

(1  -  6  A)2 

1  -  12A 

(3.14) 

/t  20  = 

V*20  = 

1 

l 

if 

1  -  1 8  A 

■ 

6 

1  -  12A 

3-stage,  scheme  of  order  5  (TG(3,5)  ) 

We  make  Z\,  Z2  and  Z3  to  be  OCA/3),  OCA/4)  and  OCA/5),  respectively,  which  leads  to  the 
foiiowing  system  of  equations: 

At  10  =  ci 
2^10  =  cj  -  2A 
6A  =  cf 
At20  +  A*21  =  ci 

2(v2o  +  ^21)  +  2fiuC\  =  Cj  -  2 A 

2  3 

6v2iC[  +  3//21  c,  =  c2  -  6Ac2 

12v2,cf  +4/i21  cJ  =  C2-  12Ac*  (3.15) 

At  30  +  A^32  =  1 

2( V30  +  K5,  +  Vn  +  //32)  =  1-  2A 

6(V31  +^32C2)+ 3//32  C2  =  1- 6A 

12(V31  +  ^2^)  +4/232  Cj  =  1  -  12A 
20(^31  +  v32 c\  )  +  5/i32  Cj  =  1  -  20A 
At2i  =0 

For  stability  purpose  we  require  additionaly  that 

V32(V2lVio-  V2oA)  =  0  (3.16) 

A  solution  to  (3.15)  and  (3.16)  is  given  below: 

ci  =  At  l  o  =  (6A)I/2 

Vio  =  2A 
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C2 

= 

J5  <28±8> 

Ml 

= 

C'\  l  <ic2  '  2  C2  *  °2  '  c2>} 

V21 

= 

S2  [  J  4  '  3  clc2  *  A<3c2  *  2cl  °2  )  1 

MO 

= 

C2  '  M 1 

v20 

= 

-  Cj  -  A  -  V2i  -  Ml^ 

M2 

= 

(6cj  -  4 Cj  ) 

V*32 

= 

0 

H  31 

= 

0 

V*31 

= 

6  •  A  *  \  ^32C2 

MO 

= 

1  -  fj.^2 

V*30 

= 

\  -  A  -  V31  -  ^32 

3-stage  imbedded  scheme  of  orden.4  (TG(3,4)  _) 

To  allow  for  an  error  control,  we  construct  a  method  in  which  a  1 -stage  2nd  order 
scheme  is  imbedded  in  a  2-stage  3rd  order  scheme  which,  in  turn,  is  imbedded  in  a  3-stage 
4ih  order  scheme.  (We  symbolically  write  TG(1,  2)  “c”  TG(2,  3)  “c”  TG(3,  4)  ). 
Toward  this  end,  we  make  Z\,  Z2  and  Z3  to  be  0(At2),  0( Ar5)  and  0{ Ar4),  respectively, 
and  we  set  c\-  =  1.  Thus,  each  stage  is  an  approximation  to  the  solution  U  at  time  tn  + 

Af  and  the  difference  ||Z/-Zj.|[|  defines  a  relative  error  for  Z/.j.  Resulting  order 
conditions  are  listed  below: 

l*\o  =  1 
2vjo  =  1  -  2A 
^20  +  ^21  =  * 

2(^20  +  v2i)  +  2^21  =  1-  2A 
6V21  +  3^21  =  1*  6A 

^30+^32-1  (3.18) 

2(  V*30  +  '/31  +  ^32  +  H 32)  “  1-  2A 
6(  1  +  V32C2)  +  3/^32  —  1  -  6A 
1 2( v*3 1  +v32  )  +  4^/32  =  1  -  12A 
Mi  =° 
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For  stability  purpose  we  require  additionaly  that 


V21  vt0  -  V20A  =  0 

V32(V2i/itO  -  ^21^10  -  =  0  (3.19) 

A  solution  to  (3.18)  and  (3.19)  is  given  below: 


M10 

= 

Ci  =  1 

V10 

= 

5  -  ^ 

Cl 

= 

1 

V20 

= 

3  (| 

4 

V21 

— • 

3  ^  Ts 

4 

8 

1*21 

- 

i  -2A  - 

/*20 

= 

-I 

32 

= 

1 

2 

V32 

= 

0 

/*31 

= 

0 

V31 

= 

-Ti 

fJ-30 

= 

0 

1 

12 

1*  30 

= 

18 


8  '  * 


3  A  (*  W  +  A  ) 


A 

1 


3 

2 


A(--fe  +  A  ) 

I  1 
8  K 


(3.20) 


Next,  using  the  original  equations  (3.1)  -  (3.5),  we  calculate  the  time  derivatives  in 
terms  of  spatial  derivatives  as  follows: 


U  t  =  -  AV 


■  Fk(Vu)k  \ 

K  ~  v  > 


(3.21) 
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(3.22) 


U,„  -  -  (AU  )_, 

rtR*'(Fu)v,M 


dF^ 

where  J?w  =  are  the  Jacobian  matrices  corresponding  to  fluxes  Fk.  In  deriving  (3.22) 
we  used  the  following  relations: 

%  =  -  (Fk{  Vu)  k  )  t 


=  -  (Fk(  Vu)  t  )k 

= 

=  -(«"(  Vu)vj\k 


(3.23) 


Also,  it  is  important  to  notice  that  UJt  can  be  expressed  in  terms  of  spatial  derivatives  of  u 
and  v  of  the  order  at  most  2,  which  can  be  effectively  handeled  by  C°  continuous  finite 
elements. 

Next,  replacing  the  time  derivatives  in  (3.7)  by  formulas  (3.21)  and  (3.22),  and 


denoting  Z.  = 


,  u 


and  Un  = 


,  we  arrive  at  the  following  system  of  equations 


Vi  -  Ov*'  1  -r\Vuj)^  }  + 

=  vn  +  KftAt  [-f*(  Vun)  k  J  +  KaAt2[-(  Rkl(Vun)  vnj)  k] 

(3.24) 

«/  ~  ^=i  (ViA*  ]  +  VijAt2[-Fk(Vuj)  k  ]  ) 

=  un  +  Ki2At  [vn  ]  +  K^At2[-Fk(Vun  )k  ] 

i  =  1,2,  ....  s 


where  the  indices  i  andy  are  used  to  denote  a  particular  stage,  the  indices  k  and  /  refere  to 
the  axis  of  a  Cartesian  coordinate  system,  comma  denotes  partial  differentiation,  and  V 
summation  convention  for  k  and  /  holds. 
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Finally,  multiplying  (3.13)  by  vector-vaiued  test  function  W  =  .integrating 

J 

over  Q.  and  integrating  by  pans,  we  arrive  at  the  variational  formulation  of  the  form: 

Find  Vj  and  ui  ,i=  1,  2, ....  s,  such  that 

J  >%  vi  dx 

~  Mi/Ar  [  J  w*  k  Fk(  Vuj )  dx  -  wTy  Fk(  Vuj  )nk  ds  ] 

-  VijAt2  [  |  wTvJc  (  Ru{  Vuj )  Vjj  )  dx  -  J  wTv  (  Ru(  VUj )  vjj  )nk  ds  ] 

£2  dil 

=  J  yvT  vn  dx 

-  Ki2At  [  I  w]  ,Fk{Vun)dx  -  I  wT  Fk(  Vun)nk  ds  ] 

•'q  v*  JXl 

-  KftAi2  [  [  wT  .  (  Rkl(  Vun)  vn  i  )  dx  -f  wT  (  Ru(Vun)  vn  i)nk  ds  ] 

JQ  ’  JdCl 

J  wl  «i  dx 

'  X/=i  I  wlvi  dx 

~  y,M  Vi jAt2  [  f  >V T  k  Fk(  Vuj  )dx  -  J  wT  Fk(  Vuj  )nk  ds  ] 

=  I  wT  un  dx 

-  KftA t  J  wTvndx 

-  KfiAt2  [  J  w]  .Fk(Vun)dx  -f  w] ’  Fk(Vun)nkds  ] 

n  '  J«9n 

for  all  admissible  test  functions  and  h»u  . 

where  n  -  ( nk )  is  the  unit  outward  normal.  Weak  form  (3.25)  is  the  basis  for  FE 
approximations. 


(3.25) 
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3.2  Taylor-Galerkin  Schemes  for  Linear  Systems 

We  now  give  an  alternative  derivation  of  TG  schemes  which  is  valid  for  linear 
systems  of  conservation  laws  with  constant  coefficients  addmiting  spectral  decomposition 
of  underlying  spatial  operator. 

Given  the  solution  Un  =  V(tn)  at  time  tn  =  nAt,  we  seek  the  next  time  step  solution 
yn+i  _  +  in  the  following  form 

[(/  -  XAp-cP-fdt2) ...  (/  -  XAt2^ldt2)]un^ 

1 - s  times - 1 

=  [/  +  *,(A)A td/dt  +  *2(,l)Ar2<?W  +  ...  +  ]un  +  0( At™+1)  (3.26) 

where 

Xjtt)  e  R,  j  =  1,  2,  ....  2s 
XfX)  =  0,  j  =  m  +  1,  wt  +  2, ....  2j 

5  =  number  of  stages 

m  <  2r  is  the  order  of  the  highest  derivative  at  the  right-hand-side  of  (3.26). 

Coefficients  Xj(X)  are  to  be  chosen  so  as  to  obtain  the  highest  possible  order  of  accuracy, 
subject  to  stability  or  other  constraints.  As  previously,  a  free  parameter  X  is  to  be  chosen 
from  stability  considerations. 

The  coefficients  ^(A)  are  determined  by  expanding  Un+1  at  the  left-hand  side  of  (3.26) 
in  a  Taylor  series  about  Un  and  equating  coefficients  of  powers  of  At  to  zero.  Accordingly, 
identity  (3.27) 

[(/  -  XAt2£ldt2) ...  ( I  -  XAt2Pldt2)][l  +A td/dt  +  ^Af2<?W  +  ...  +  ■~At2s^s/dt2s  ]un 
1 - 5  times - 1 

=  [i^X\KX)AM  +^2(A)A t2cP-lch2  +  ...  +  j£25(A)A/1?<^-7<3f25  ]f/n  +  <9(A/™+1) 

(3.27) 


leads  to  the  following  relations  for^(A) 


*i(A)  - 1 

xm  =  £  -  (0  A 

*3W-  5?  -  (0A 

z,W  =  jr  -  5t(')  a  +  (*)  A2  (3.28) 

««-s-i(0a  +  ii(0A2-G)aI 

aw-MC)***©*-  (O*1 

etc. 

where  =  -,J!^  is  the  Newton  binominal  symbol. 

Constructed  in  this  way,  an  5-stage  scheme  is  of  the  order  m  <  2s.  Furthermore,  for 
particular  values  of  X  it  will  be  of  the  order  m+1,  and,  in  fact,  this  is  the  highest  attainable 
order  of  accuracy.  Typically  we  choose  m  -  2s,  i.e.,  an  5-stage  scheme  of  the  order  2s  (or 
2s  +  1).  On  the  other  hand,  the  choice  m  <  2s,  when  some  of  the  coefficients  X]  are  a 
priori  specified  to  be  zero,  may  be  advantageous  from  the  stability  point  of  view,  e.g.,  by 
setting  X2 s  =  Xzs-\  -  0»  the  resulting  5-stage  scheme  will  be  of  the  order  25-2  (or,  at  most 
25-  1). 

Remark:  the  scheme  with  m  <  2s  ,  e.g.,  m  =  2s  -  1,  where  %2s  s  should  not  be 
confused  with  the  scheme  with  m  =2s  and  X  being  so  chosen  that  %2s  ~  The  former  is 
of  the  order  2s  -  l  while  the  letter  is  of  the  order  25.  It  is  only  for  this  particular  value  of  X 
that  *h  schemes  are  identical.  H 

EXAMPLE  3.1 

The  3-stage  scheme  of  the  order  6  is  of  the  form: 

[/  -  XAt26>1/dt2Yun+*  =  [/ 

+  At  dldt 

+  (  \  -  3A  )  At2  P/di2 
+  (  -  3A  )  At3  aw 
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+  -  §A  +  3A2)Ar4<?W 

+  (f4o  '  2  A  +  3A2)  Ar5  oW 

+  (  '  8  A  +  2  A2  '  ^  A/*  ]^"  (3.29) 

with  the  local  truncation  error  Et  satisfies 

[/  -  AAr2<?W]3£,  =  [( 5^40  -  ^  ^  +  i  A2  -  A3)  Ar7  d7ldt7)uH+0( Ar8)  (3.30) 

i _ <?,(A) _ i 


Setting  A  -  1.41218087134444  yields  %6  =  0  but  the  resulting  scheme  is  still  of  the  order  6. 
By  choosing  A  to  be  a  zero  of  e,(A)  (e.g.,  A  =  0.444797521031781),  formula  (3.29)  has 
order  7. 

The  3-stage  scheme  of  the  4th  order  is  of  the  form: 

[/  -  AAr2^2/^2]V/!+1  =  [/ 

+  A tdidl 

+  ( ^  -  3A )  A/2  cP-ldt2 
+  ( -  -  3A )  Ar3  cP/dP 

+  ( 24  •  |a  +  3A2)  (3.31) 


In  order  to  implement  formula  (3.27)  efficiently,  it  is  necessary  to  factorize  the  right- 
hand-side  operator  into  the  product  of  “quadratics”  as  follows: 


(/  +  £,(A)A tdldt  +  ;£2(A)Ar2c£/cfr2  +  ...  +  ^(A)A ficps/dt25  ) 

=  (/  +  /r,(A)Ar<7/e?/  +  v,(A)A p-cP-ldt2)  ...  (/  +  fi%(X)Atd/ck  +  v^(A)A t2cP/dt2) 
1 - — - - — — — — s  times - 1 


(3.32) 


Indemity  (3.32)  leads  to  the  following  system  of  nonlinear  algebraic  equations  for //y  and  , 
7=  1.2 . s: 
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Mi+M2  +  -  +  =  1 

vt  +  v2  +  ...  +  vs  +  /i^2  +  ...  +  MslM,  =  *2<A) 


viv2  -v,  =  *2,<A) 


(3.33) 


Factorization  (3.33)  exists  and  is  unique.  The  coefficients  of  the  factorization  can  be  found 
numerically  for  a  given  value  of  parameter  A  and  are  given  in  an  Appendix. 

Next,  we  use  relations  (3.21)  and  (3.22)  to  express  time  derivatives  in  term  of  spatial 
derivatives.  For  linear  systems  with  constant  coefficients,  these  equations  take  a  particulary 
simple  form: 


V  t  =  -  AU 

V  M  =  A2U  (3.34) 

Replacing  the  time  derivatives  in  the  left-hand  side  of  (3.27)  and  the  right-hand  side  of 
(3.29)  by  formulas  (3.34)  and  introducing  auxiliary  stages  Y( ,  we  arrive  at  the  following 

system  of  equations 

K,  -  AAf2A2K1  =  Un- fi]AtAUn +  \\At2A2Un 

Y2-XAt2A2Y2  =  Yl-u2MAYl  +  v2At2A2Y  x 

s  times  (3.35) 


[_  t/n+1  -  XAt2A2Un+i  =  -  nsAt AYs  X  +  vsAt2A2Ys_l 

In  deriving  (3.35)  we  used  the  fact  that  the  operators  [ I  -  AAr2A2]"1  and  [/  -  ^Ar  A  + 

v At2 A  2  ]  commute.  (In  general,  [f{A  )og(/t)](f/)=  V  00  f((0  )g((0  )dP  (U)  = 

[£(A)o/t4)](f/),  where  f{A)  is  a  bounded  operator,  con  are  the  eigenvalues  of  A  and  dPn 

are  the  associated  eigenprojections;  see  also  Section  2). 

As  can  be  seen.  (3.35)  is  a  sequence  of  s  linear  elliptic-like  PDF’s  with  the  same  left- 
hand  side  operator.  Moreover,  since  operator  A2  has  zero  off-diagonal  terms  : 
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<4 2  = 


-Fk(m) 


-Fk(V(-)) 


(3.36) 


a  typical  “one  stage  problem”  is  splitted  into  v-equations  and  u-equations,  which  can  be 
solved  independently  and  in  parallel  (the  presence  of  A  at  the  left-hand  side  of  (3.35) 
would  have  destroyed  this  property).  In  addition,  the  diagonal  terms  of  A 2  are  identical 
which  means  that  the  operators  defining  the  left-hand  side  of  each  stage  of  (3.35)  are 
identical  for  both  v-  and  iz-equations  (and  for  all  stages). 

Finally,  multiplying  (3.35)  by  vector-valued  test  function  W,  integrating  over  Q  and 
integrating  by  parts,  and  denoting  F0  s  Un  and  Y s  =  Un+l  ,  we  arrive  at  the  variational 
formulation  of  the  form: 


Given  YQ  e  X 


:ind  Yt  =  je  X  ,  i  =  1,  2, ....  s, 


such  that 


a(H’v,v/)  -  XA t2[b(wv,v;)  -  ,  vt)J 

=  a(wv  ,  v(.j)  +  /ri-Ar[6(H’v  ,  u,-.,)  -  b^Wy  ,  uM)J 

+  v.  A t2[b(wv  ,  vi _,)  -  b^Wy  ,  V..,)] 

a(wu,  ut)  -  AAr2[6(M-tl ,  u-)  -  b^Wu,  «.)] 

=  a(wu  ,uiA)  -  /iiAtfl(HP-lvM) 

+  v.Apibiwu  ,  u(.  j)  -  bf{wu,uiX)) 


(3.37) 


for  all  admissible  test  functions  W  = 


-  c:  )■ 


Here  : 


a(w,  u) 

|  wru  dx 

b(w,u)  ~ 

w  T.  Rklu  ,dx 
Jn  * 

m 

bjiw,  u)  = 

wTn.Rktuds 

(3.38) 
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jpk 

Jacobian  fluxes  Ru  =  are  now  constant  matrices,  and  nk  is  the  unit  outward  normal. 

Replacing  X  in  (3.37)  with  a  finite  dimensional  subspace  Xh  of  X  ,  we  arrive  at  the 
fully  discretized  problem. 


Given  rh>0e  Xh 

Is  ■i  =  1-2 . *• 

such  that 

a{yvy,vhi)  - 

AAr2[&(w„  ,  vhi)  - 

bfiWy  ,  %  ■)] 

=  a(wv  ,  vw.j)  + 

+  V;  At2[b(wv  ,  %•.,)  -  b^wy  ,  vhM)] 

(3.39) 

a(wu,uhi)  ~  XAt2[b(wu  ,uhi)  -  b^Wu  »«/,,,)] 

=  a(»u  .  uh  iA)  -  a(wu  ,  vhi  i) 

+  v,.  A t2[b(wu  ,  uh  i_ j)  -  bfiwu  ,  «w.,)] 

(w  v  'l 

for  all  admissible  test  functions  W  =  e  X. 

J 


Formula  (3.39)  is  referred  to  as  an  5-stage,  high-order  Taylor-Galerkin  (TG)  method. 


4  A  Linear  Stability  Analysis 


We  shall  now  restrict  ouerselves  to  the  special  case  of  the  eqautions  of  linear 
elastodynamics.  Thus,  we  consider  the  problem  of  the  form 


d^u  tm 

w  +  a“  * 

du 

U  =  uo  >df 


0 


=  v , 


t  >  0 
t  =  0 


(4.1) 
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where  the  operator  8  is  defined  by  eqns  (2.10)  -  (2.16).  It  is  convenient  to  convert  (4.1) 
into  the  1st  order  system  of  the  following  form 


where 


{U  r  +  iA  U  =  0  t>  0 

U  -  U  q  r= 0 


dcf  f  0  8  \ 

AU  =  -  i 

.  -I  0 


u 


(4.2) 


(4.3) 


U  is  a  group  variable,  U 


^  j  ,  UQ  specifies  initial  conditions,  U0 


,  and  i  is 


the  imaginary  unit.  The  form  of  (4.3)  is  similar  to  (3.4)  except  the  multiplier  i,  which  is 
introduced  to  make  the  operator  A  self-adjoint. 


The  multi-stage  Taylor-Galerkin  method  for  solving  (4.2)  reduces  to  a  sequence  of  5 
linear  variational  boundary-value  problems  of  the  form 


Given  Un  e  X 

Find  Kj,  F2,  ...,  F^,  Un+i  e  X  such  that 

A(W,  F,)  +  AA t2B(W,  Ft) 

=  A(W,  Un)  +  /i,A tC(W,  Un )  -  v,A t2B(W,  Un) 

A(W,  Y2)  +  AA t2B(W,  Y2) 

=  A(W,  F,)  +  /t2A/C(W,  F,)  -  v2At2B(W,  F,) 

s  times  (4.4) 

A(W,  Un+]  )  +  XAt2B(W ,  Un+l  ) 

-  A(W,  F (  , )  +  fisA tC(W,  Fs  |)  -  vsA t2B(W,  FM) 

for  all  admissible  test  functions  We  X 
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where  X  =  Xv  x  Xu  —  H  x  D{ SI),  the  bilinear  (sesqulinear)  forms  A,  B  and  C  are  defined 
by 


A,B,C  :  (Xvx  Xu)x(XvxXu) — >  C, 

A(W,U)  =  (V>W)HxH 

B(W,U)  =  ( AU,AW)HxH  (4.5) 

C(W,U )  =  -i(U,AW)HxH 

and  the  inner  product  (•,  -)H  x  H  is  defined  by 

{V,W)HxH  =  ([v,  u],  [Wy,  WU])H  x  H 

(4.6) 

=  (v,W*)H+(u,Wu)H 


Each  stage  of  (4.4)  defines  a  linear  operator  Tj  from  X^=Xn  Rg(Tj )  into  itself; 


7> :  - 

zrTizj-\ 


(4.7) 


and  the  composition 


T  =  7>7V,o  ...  or, 


(4.8) 


defines  a  transient  operator  T  taking  an  approximate  solution  Un  at  time  level  f„into  Un+l 


T  :  X , 


Un+1  =  TUn 


(4.9) 


It  should  be  noted  that  bilinear  forms  forms  A  and  B  are  symmetric  (Hermitian) 


A{W ,  U)  =  A(U,  W) 


5(W,  V)  ~  B(U,  W) 


(4.10) 
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and  bilinear  form  C  is  skew-symmetric  (skew-Hcimitian) 


C(W,U)  =  -  i(U,AW)HxH 
=  -  i(AU,  W)H  x  H 

=  -  i(U,A\V)  HxH  (4.11) 


-  i(U,  A  IV)  H  x  H 
=  -  C(U,  W) 

Moreover,  the  bilinear  form  appearing  at  the  left-hand  side  of  (4.4)  satisfies  the  inf-sup 
stability  condition 


sup  U(W,  t/)  +  AAr2fi(W,  f/)  |  £  {U,  U)HxH  +  \M2(AU,AU)HxH  (4.12) 

and  defines  a  natural,  energy  norm  for  problem  (4.4).  Condition  (4.12)  implies  also  that 
the  variational  problem  (4.4)  posesses  a  unique  solution  Un+l  and  guarantees  the  usual 
convergence  properties  for  variational  boundary-value  problems. 

To  investigate  stability  properties  of  (4.4),  we  need  to  estimate  the  eigenvalues  of 
transient  operator  T,  or,  equivalently,  the  eigenvalues  of  the  the  following  eigenprohlem 


Find  an  eigenpair  (f2,  0  *  U  e  X  ) 
and  Tj,  Y2,  ...,  Y  s Ae  X  such  that 

A(W ,  Fj)  +  AAr2fl(VV,  V,) 

=  A(W,  U)  +  /i,ArC(VV,  U)  -  V^BIW,  U) 

A(W ,  Y2)  +  XAt2B(W,  F2) 
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=  A(W,  F,)  +  p2A tC(W,  Fj)  -  v2A/2fi(W',  F,) 


(4.13) 


^  times 


fl[A(W,C/)  +  AAf2B(VF,  (/)] 

=  A(W,  YsA)  +  ^A/C(W\  F,.,)  -  vsAi2B(W,  Ys_ j) 


for  all  admissible  test  functions  VFe 


Toward  this  end  we  notice  that  the  underlying  operator  A  is  self-adjoint  (in  the  complex 
sense)  and  has  a  pure  point  spectrum 


o(A)  =  { 0,  ,  (O, ,  ,  (02  ,  2  , ... } 


0  <  G),  <  GJ2  <  ...  <  con  — >  oo  ,  0)n  e  R+ 


with  corresponding  orthonormal  eigenfunctions  [U  }  such  that 


(4.14) 


T-l°° 

*  =  a*dp« 
'  -  dP0  +  X°n-  JP" 


(4.15) 


(4.16) 


Un)  =  5mn 


(4.17) 


where  dPn{‘)  -  A(-,  Un)Un  ,  dP0  is  the  projector  onto  N(A),  co  n  =  ~(on  and  U  n  =  -  U n 
By  taking 


u  =  dp0(v)  +  y~  dpnm 


(4.18) 


zi  - dp^) - 


(4.19) 
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W  =  un  €  N(A) 


(4.20) 


eigenproblem  (4.13)  can  be  reduced  to  the  following  form: 


Find  an  eigenpair  (£2n  ,  0  ±U  e  A'  ) 
and  Yy  Y2 , j€  J9(A)  such  that 

A(Un,Y{)  +  ZAt*a£MUH,Yl) 

=  A(Un ,  U)  -  i^AtconA(Un,U)  -  v,A t2co2nA(Un  ,  U) 

A(Un,Y2)  +  XAt2co2nA(Un,Y2) 

=  AiUH,Yx)-ifi7Ata)nA(UntYJ  -  v2At20)2n  A(U  n ,  K,)  (4.21) 

j  times 

[/»(£/„ .  t/)  t  AAA<>*,MJ/„  ,  U)\ 

=  -W„  •  I'm)  -  HtA>o>.MV, .  I",,)  -  vsAr^  -4(t/„ ,  V,.,) 


for  Un<=  X 


From  (4.21),  we  can  obtain  Qn  , 


Q 


■.  -  n,i, 


1  -  v>(A)Ar2^  -  //i>(A)Ar 
1  +  AAr2^ 


(4.22) 


and  its  modulus 


ifl  p  =  n5  0  - 

"  A1'-‘  ( 1  +  AArto*  )2  ' 

By  choosing  U  -  VV  e  A 1(A)  it  is  easily  seen  that  the  corresponding  eigenvalue  of  T, 
is  equal  to  unity: 

420  =  1  (4.24) 
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Similarly,  the  eigenvalues  of  the y-th  stage  eigenvalue  problem  : 


Find  an  eigenpair  {Qjn,  0  Z  e  X  )  such  that 

DJn[A(W,Z)  +  Ut2B(W,  Z)\  (4.25) 

=  A(W,  Z)  +  HjAiC(W,  Z)  -  VjAt2B(W,  Z ) 

for  all  admissible  test  functions  We  X 


read 

_  1  -  v,(^At2o)2n  -  iH,(X)At  an 
*'n  1  +  A  A  r2<y2 

(4.26) 


and,  therefore  (4.22)  and  (4.24)  can  be  written  as 

a.  =  fl;=,  nj,  «=«■'.  •••  (4.27) 

In  other  words,  the  eigenvalues  of  the  transient  operator  T  are  products  of  the  eigenvalues 
of  the  component  operators  T '■ . 

To  show  unconditional  stability  of  TG  methods,  we  need  to  show  existence  of  a  A  such 
that 


\Q^  <1,  n  =  1,  2,  ... 

V  Ate  R+ 


(4.28) 


Thus,  in  order  to  find  a  A  producing  unconditionally  stable  Taylor-Galerkin  method,  we 
need  to  solve  the  following  problem: 

J~ Find  A  €. it+  such  that 

(1  -  v-(A)A t2(o2n)2  +  (fij(h)At  con  )2 

(1  +  AA t2(o2n  )2  ~  1 
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Mi  +M2  +  -+Ms  =  1 

v,  +  v2  +  ...  +  v5  +  hxM2  +  +  Ms.iMs  =  c2(X) 

!  (4.29) 

VlV2  ...V,  =  C2sa ) 

for  all  At2®*  e  R+ 

Unfortunately,  (4.29)  can  be  solved  analytically  for  s  =  1  only.  In  this  case,  it  reduces  to 
finding  X  such  that 


(l  -  X)At2(o2n)2  +  At2(o2n 

( 1  +  XA t2col  )  2 


<  1 


(4.30) 


which  is  satisfied  for  all  At  provided  X  >  ^ .  (This  result  was  proven  in  an  alternative  way 
by  Demkowicz  et  al.  in  [1].)  For  s  >  1  we  compute  and  v-  for  a  given  X  and  verify 
stability  criterion  (4.30).  The  results  of  numerical  computation  of  the  ranges  of  X  which 
satisfy  (4.29)  are  summarized  in  Tables  4.1  and  4.2.  In  addition,  Figs.  1-7  show  the 
variation  of  \QJ  as  a  function  of  A r2<y2  for  X's  lying  at  the  boundary  of  stability  regions 

listed  in  Tables  4.1  and  4.2. 


Table  4.1 


s 

m 

X 

1 

2 

[4  *  °°) 

2 

4 

[0.47048,  00) 

3 

6 

[0.695,  00) 

4 

8 

[0.91943,  00) 

5 

10 

[1.15,  00) 
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Table  4.2 


5 

m 

A 

3 

4 

[0.22,  oo) 

4 

6 

[0.29,  oo) 

5 

8 

[0.38,  oo) 

It  is  interesting  to  check  whether  the  5-stage  methods  of  order  2s+l ,  i.e.,  the  methods  of 
the  highest  attainable  order  of  accuracy,  are  unconditionally  stable.  Table  4.3  gives  the 
largest  values  of  A  producing  5-stage  methods  of  the  order  25+1.  Unfortunately,  all  those 
A’s  are  not  in  the  range  of  unconditional  stability,  and,  therefore,  the  5-stage  (25+l)th  order 
methods  can  be  at  most  conditionally  stable. 


Table  4.3 


5 

m 

A 

1 

2 

l 

6 

2 

4 

1  +  1.J2L 

6  2  V  270 

3 

6 

0.444797521031781 

4 

8 

0.583260771000355 

5 

10 

0.721627161040466 

Next  we  check  stability  properties  of  the  5-stage  25-th  order  methods  with  A  so  chosen  that 
coefficient  c2s  =  0.  Table  4.4  lists  the  smallest  of  such  A’s  and  Figs.  8-11  show  the 
corresponding  variations  of  eigenvalues  as  a  function  of  Ar2ct)^  .  As  expected,  all 
those  A’s  produce  unconditionally  stable  Taylor-Galerkin  methods. 

Table  4.4 

5  m  A 


2  4  0.956435464587639 
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3  6  1.41218087134444 

4  8  1.86773692752801 

5  10  2.32321418507049 


Moreover,  this  class  of  methods  has  the  property  that  the  eigenvalues  Qn  tend  to  zero  as  At 
goes  to  infinity 

lim  1^1=0  (4.31) 

At  — >oo 


This  is  easily  seen,  since  from  (4.23)  we  have 


lim 
At  — >  oo 


l"J 


|v,v2...vj 


(4.32) 


and  one  of  the  Vj  is  equal  to  zero.  In  fact,  all  the  methods  with  m  <  2s  will  have  property 
(4.32). 

Equation  (4.32)  suggests  also  a  criterion  of  selection  of  A  for  methods  with  m  =  2s. 
Namely,  of  all  the  A’s  guaranteeing  unconditional  stability  we  choose  the  one  for  which 
|vt  v2  ...vs  |/A*  is  the  largest,  i.e.,  we  choose  the  A  which  will  produce  the  method  with 
least  possible  numerical  dissipation.  These  optimal  values  of  A  as  well  as  corresponding 
asymptotic  values  of  ]Qn\  are  given  in  Table  4.5.  As  can  be  seen,  all  the  methods  but  s  =  1 
are  dissipative.  Note  also  that  the  optimal  A’s  lie  at  the  boundary  of  stability  regions  listed 
in  Table  4.1. 


Table  4.5 


s 

m 

A 

l^i  v2  ...vJ/A* 

1 

2 

l 

4 

1 

2 

4 

0.47048 

0.937 

3 

6 

0.695 

0.904 

4 

8 

0.91943 

0.887 

5 

10 

1.15 

0.868 
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Figure  12  illustrates  the  decomposition  of  the  eigenvalue  of  the  transient  operator  T  into 
the  product  of  the  eigenvalues  of  the  component  operators  Tj  for  the  4-stage  8th  order  TG 
method  (compare  eqn  (4.27)).  As  can  be  seen,  the  moduli  of  eigenvalues  of  Tj  are  bounded 
but  they  are  not  necessarilly  less  than  unity  for  all  the  component  operators.  (This  property 
also  holds  for  other  unconditionally  stable  schemes).  In  general  we  have 

\QjJZCj  ,  7=1,2 .  5  (4.33) 

where 

Cj  =  ess  sup  |  (4.34) 

0<A/26)^<oo 


usually  Cj  =  max  { 1 ,  |  v  |/A  } . 


Thusfar  we  investigated  stability  properties  of  TG  methods  for  Xx  being  infinite¬ 
dimensional  and  we  showed  the  existence  of  ranges  of  the  stability  parameter  A  for  which 
the  spectral  radius  of  T,  r(T),  is  less  than  or  eqal  to  unity.  To  show  unconditional  stability 
of  TG  methods  for  finite  dimensions,  i.e.,  when  T  is  restricted  to  Xh  (Xh  c:  X j,  dim  Xh 

<  oo),  we  will  show  that  the  spectral  radius  of  T  satisfies  the  following  maximum  principle: 


where 


sup 
C/e  AT, 
U*  0 


BX{TV  ,  TV ) 
BjU  ,  U) 


=  r2(T) 


B{(U,  U)  =  A(U,  U )  +  AA t2B(U,  U) 


(4.35) 


(4.36) 


is  the  bilinear  form  appearing  at  the  left-hand  side  of  (4.4),  and 


r(T) 


-  sup 


l«J 


(4.37) 
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Indeed, 


B^U.U)  = 

=  B, (<//■„(£/)  +  dP„(V)  , dP0(U)  +  ) 

-fl,(dP0(t/)  +  '^4Z.^U’U^fJr,  -dP0(U)  *  ’ZZ.„A(-lJ-U^Vn) 

=  A(dP 0(t/),  dP q(U)  )  ^  -  zr,„  ) 

+ aa,2/'(£2.„  “v^-  ) 

~A(dP0(U)  ,  <tf>0({/)  )  +  ^  |A(t/.  t/„)|2(l  +  AAl2|<aJ2  )  (4.38) 

i 

b^ti/,  rt/)  = 

=  Bt(r  (dP0(U)  +  ^i^dPJfl)),T(dP^V)  +  «l/))  ) 

=  fl,(  <iP0(l/)  +  ,dPj.U)+YJZ.J2^V-V»'>V»  ) 

-  4(d/>0(£7)  ,  <//>„({/)  )  +  *(£“  U")U"  ’  Un>Un  ) 

+ aa'2/i(X2-„  «.)w. ) 
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=  A(dP0(V)  ,  dP0(U)  )  +  ]T~  ^  | Qn\2  \A{U ,  t/„)l2(l  +  Mt2|<uJ2  ) 

<  A(dP0(U)  ,  ^0(^)  )  +  supl^J2  ^  M*7’  f/")|2(1  +  AA/2I"J2  ) 

<  max {  1,  sup|l2j2}  \^A{dP0{V)  ,  dP0(U)  )  +  ^  \A(V,  Vn)\2(\  + 

=  max  (|fl0|2,  suplflj2)  U) 

=  r2(T)  BJU,  V) 

Thus, 


sup 
Uzxx 
U*  0 


BX(TU  ,  TV) 
~BX(V  ,  V) 


<  ^(T) 


Conversely, 


l«J2  = 


Bx(TVn,TVn) 


BATV , TV) 
sup  — L- — - — - - - 

UeX,  BX(V,V) 
U*  0 


from  which  (4.35)  follows. 

Consider  now  the  finite  dimensional  counterpart  of  eigenvalue  problem  (4.13): 

Find  an  eigenpair  {Qh%  0  *■  V h  e  IfJ 
and  7,h  x,Zh  2 . Zh  f  ,e  Xh  such  that 

A(\V,  Zh  ,)  +  AAr25(VV,  Zh  ,) 

=  A(\V,  Vh)  +  /i,A/C(lV,  Vh)  -  vxAt2B(\V,  V h) 


Mf2|«J2  )] 


(4.39) 


(4.40) 


(4.41) 
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(4.42) 


A(W,  Zh  2)  +  AA t2B(W,  Zhi2) 

=  A(W,  Zh  ,)  +  ft2AtC(W,  Zh  ,)  -  v2At2B(W ,  ZK  ,) 

s  times 

nh  [aiw,  uh)  +  xm2B(  w,  uh)] 

=  A{W,Zh  ^)*lx^,C(W,Zh,s.O  -  vi&i1B(W,Zh  l  ]) 
for  all  admissible  test  functions  W  e  Xh 


Since 


max|i2AP  < 


sup 
UeXh 
U  *0 


B  }(T  U  ,  TU) 
BX(U  ,  U) 


UeXx 
U  *0 


B  X(T U  ,  TU) 
B  j  (U  ,  U) 


(4.43) 


=  r\T) 


unconditional  stability  of  the  TG  scheme  for  infinite  dimensional  levels  (i.e.,  r2(T)  <  1  V 
At  e  R+)  implies  its  unconditional  stability  at  finite  dimensional  levels,  independently  of 
the  approximation  used  with  respect  to  space  variables  (any  finite  element  mesh). 


5.  Adaptivity  (Linear  Problems) 


5.1  Error  Estimation  for  High-Order  TG  Methods 


As  a  starting  point,  we  consider  the  estimation  of  the  spatial  approximation  error 
commited  at  a  typical  stage.  The  y'-th  stage  problem  reads 
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Given  Z.  .  ,€  X.  „  find  Z,  „  -e  Xh  „  such  that 

n.p.j’l  K  p  n,p,j  n,  p 

biW>zh.p:?  -  VW*Z*.p;y-i>  C-D 

for  all  admissible  test  functions  We  X. 

n,  p 


where 


fi,(W,  U)  =  A(W,  l/)  +  AAr2S(W,  t/) 

L;(W,  f/)  =  A(W ,  */)  +  ^ArC(W,  U)  -  vAr2fl(W,  l/) 


(5.2) 


X h  c  A’,,  dim  Xh  p<  oo  and  indices  6  and  p  refer  to  the  use  of  arbitrary  h-p  finite 
element  meshes,  with  locally  varying  mesh  size  h  and  spectral  order p. 

Assuming  that  there  is  no  error  in  Zh  ,  we  consider  the  enriched  space  Xh  j 
corresponding  to  the  same  mesh  but  with  local  order  of  approximation  uniformly  increased 
by  one  and  we  define  the  relative  error  as 

ej"  =  z».p*w  -  ZKp.i  <5.3) 


where  Z^p+l.j  is  the  enriched  space  solution. 

We  estimate  ||^l)||£  using  the  element  residual  method  in  the  form  [5] 


sup 
We  X 


I B 


<glZ».,:)>_ 


h,  p+\ 


\m. 


=  \\Z 


h.p+l;j 


-  Z 


K  p 


>lli 


(5.4) 
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where 


*  II* Ilf-  is  the  global  energy  norm  defined  as 

l!tf|I£  =  BX(U,  U)v2  (5.5) 

•  TL.  •  is  the  element  error  indicator  function  evaluated  as 

Ki^'n  <5-6> 

with  Bx  ^.the  element  contribution  to  the  (global)  bilinear  form  Bx  and  the  <pK  j  the 
element  indicator  function  which  is  the  solution  to  the  local  problem 

Find  <pK  ■  €  X®  p+l(K)  such  that 

(5.7) 

=  Rj  W  v  We  <p+1(K) 

where  Rj  is  an  appropriate  residual  corresponding  to  element  K  and  JCh  j(K)  is  the 
kernel  of  the  h-p  interpolation  operator  defined  on  the  element  enriched  space 
*H.P+ 1(*0- 

For  all  details  we  refer  to  [5]. 

Next  we  consider  the  error  Ej2)  which  is  caused  by  the  error  in  the  previous  stage 
solution  Zh  j and  enters  (5.1)  through  the  right-hand  side.  Formally  Ej  is  computed 
as 

E?  -  TH.p.UjE,A  (5.8) 

where  Th  p+l.j  is  the  operator  taking  the  previous  stage  solution  Zh  ,  into  Zh  j  . 


T h, p+l:)  '  ^ h. p+ 1 


L h,  p+l 


7  =  T  7 

h,p+\.j  ,h,p* \\j  p+ 1 ; j- 1 


(5.9) 
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which  can  be  proved  in  exactly  the  same  way  as  the  maximum  principle  (4.35). 

Thus,  the  spatial  approximation  error  accumulated  up  to  the  y'-th  stage  is  given  by  the 
following  recurrence  relation 
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IIE;II£  s  IIe!°II£  +  IIEfll£ 

(5.13) 

S  IIE,")II£  +  q  ||e,,||£ 

By  using  (5. 13}  at  each  stage,  we  can  estimate  the  spatial  approximation  error  commited  at 
a  typical  time  step  U1} - »  Un+]  .  Accordingly, 

P  P 

IIE,II£  *  IIE^’IU  +  cJej.,||£ 

I|E,.iII£  s  ||£''»||£  +  C„,  l|Ej!£  (5.14) 

IIe,II£  -  IIe'/’IU 

and,  therefore, 

iiOi* H  ifci u 

^  cs  Ki  (cs-2  ...  q  ll£(,l)ll£  +ll£2)ll£  )  +  !I4i>!'£  )  +  -  +  ll^!|llc )  + 

(5.15) 


where  it  is  assumed  that  there  is  no  error  in  the  previous  step  solution  Vnh  .  It  can  also  be 
shown  that  |j£^  ||£  is  a  bound  for  the  relative  spatial  approximation  error  commited  at  a 
typical  time  step: 


ll<ll£  *  || V 


n+1 
h,p+ 1 


(5.16) 
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respectively. 


*h.p+l  • 


To  account  for  the  total  error  commited  at  a  typical  time  step,  it  is  also  necessary  to 
consider  the  temporal  approximation  error.  We  estimate  the  relative  temporal  approximation 
error  as  follows: 


K\\e  =  \Kl  ■  <5J8> 

where  Un*p  is  the  solution  obtained  by  using  an  s-stage  m-th  order  TG  scheme  which  is 
referred  to  herein  as  the  basic  solution  scheme  and  ^n^p  is  the  solution  obtained  by  using 
the  (s+l)-stage  (m+2)th  order  scheme  (with  the  same  stability  parameter  A  and  the  same 
time  step  At  )  which  is  referred  to  as  the  auxiliary  solution  scheme.  Notice  that  the 
definition  of  auxiliary  scheme  is  unique.  Since  the  bilinear  form  defining  the  left-hand  side 
of  the  basic  method  is  the  same  as  that  defining  the  left-hand  side  of  the  auxiliary  method 
(since  A  and  At  are  unchanged),  the  computation  of  and,  hence,  ||£*  ||£  is 

economically  feasible. 

Thusfar  we  estimated  the  approximation  error  1|^  commited  at  a  typical  time 

step  and  consisting  of  spatial  approximation  error  and  temporal  approximation  error  : 

l|£<l,J%  -  l|E"J£  +  <5-l9> 

To  obtain  the  error  for  the  whole  evolution  process,  it  is  necessary  to  consider  the  error 

(2)  f\ 

||EW‘  ||£  which  is  introduced  by  the  error  in  the  previous  time  step  solution  Unh  . 
Formally, 

E{2U  =Thp+lEnl  (5.20) 

where  Th  p+x  is  the  transient  operator  taking  the  previous  time  step  solution  {/*  t 
into  and  EnA  is  the  error  in  t/”  .  (More  prccisly,  En  1  is  the  error  accumulated 

up  to  the  (rt-l)th  time  step). 


We  have  the  following  estimate: 


HE®'"  |||  =  || Th  p^  E"''||| 


=  Hn.^tll£ll£ 


|2 
I E 


sup 

V*XKp+ 1 

U  *  0 


B^TV,TU)  , 
BX(U,U)  11/1  1 


(5.21) 


< 


sup 
1/eX, 
V  *0 


BX(TU  ,  TU) 
BX{U  ,U) 


Ill’ll 


2 

E 


=  r2(T)  \\En% 


<  ||En  l 


2 

E 


where  it  is  assumed  that  r(T)  <  1. 

Thus,  the  approximation  error  accumulated  up  to  the  /i-th  time  step  is  given  by  the 
following  recurrence  relation 


l|ElE  <  I|e,1)'',IIe  +  l|£<2)''1  IIe 
s  ||e(,)-"|Ie  +  ||e"’||£ 


(5.22) 


Inequality  (5.22)  provides  a  basis  for  the  control  of  the  quality  of  results  for  the  whole 
evolution  problem. 


5.2  Adaptive  Strategy  and  Computational  Considerations 

We  neglect  the  temporal  approximation  error  (inasmuch  as  we  use  high-order 
schemes);  we  use  a  constant  time  step  A t  and  employ  a  simple  strategy  based  on 
equidistribution  of  spatial  approximate  error  for  each  time  step.  The  equidistribution 
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strategy  means  simply  that  a  number  of  elements  with  the  largest  error  indicators  are 
refined.  The  use  of  a  constant  At  is  permissible  because  of  the  unconditional  stability  of  TG 
schemes. 


Formally  the  algorithm  is  as  follows: 

Step  1  Read  in  data  (geometry,  initial  conditions,  etc.). 

Specify  time  step  At,  number  of  time  steps  NSTEP  and  the  target  error  for  the 
whole  evolution  problem  rjt . 

For  each  time  step  Un. - >  Un.+X  : 

p  «, p 

Step  2  Specify  the  target  error  for  “one  time  step  problem”  7]&  ,  e.g.,  =  ^s^Fp 

Step  3  Solve  for  Vn^x  ,  save  the  internal  approximations  Zh  p  ■  ,j  =  1,  2, ...  s- 1.  Use  a 
direct  solver  or  an  iterative  solver  using  Vnh  and  the  previous  time  step  internal 
approximations  as  starting  vectors. 

Error  estimation 

Step  5  For  each  element  K  determine  error  indicators  via  (5.6) 

Step  6  Determine  the  global  error  commited  at  this  time  step: 

I|E(,)’',He  -  iie;ii£ 

where  ||J^IIC  is  estimated  via  inequality  (5.15). 

Step  7  Check  the  error.  If  ||£(1),n  ||£  <  rjA/ ,  then  go  to  Step  13. 

Step  8  For  each  elemern  K  determine  error  indicators  for  mesh  refinements  t}*.  Identify 

the  largest  admissible  error  indicator  . 

Step  9  Create  a  list  of  elements  such  that  j]K  >  rj ^  . 
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Step  10  Check  the  number  of  elements  the  list.  If  it  is  too  small  (zero  in  particular), 

decrease  t]aci  and  return  to  Step  9. 

Step  1 1  Refine  elements  K  from  the  list  created  at  Step  9. 

Project  Unh  p,Zh  p.j,j  ~  1,  2,  ...  $-1,  and  t/"+1  onto  enriched  space. 

Step  12  Solve  for  U"+p  using  iterative  solver  and  coarse  mesh  solutions  Unh  p  ,  Zh  p.- 
j  =1,2, ...  f-1,  and  Un£p  as  starting  vectors.  Save  the  fine  mesh  solutions. 
Set  “fine”  =  “coarse”  and  return  to  Step  5. 

Step  13  For  each  element  L,  determine  error  indicators  for  mesh  unrefinements 
Identify  the  largest  admissible  error  indicator  . 

Step  14  Create  a  list  of  elements  L  such  that 

•  TfccO.Oltj  ^ 

•  $L<lad 

Step  15  Unrefine  all  elements  L  from  the  list  created  at  Step  14. 

Project  Unh  ,  Zh  p.]  tj  =  1,  2,  ...  $-1,  and  Un+p  onto  unenriched  space. 

Step  16  Set  Unh  =  Un^p  and  go  to  Step  2. 


The  element  indicator  for  mesh  unrefinements  %L  is  defined  to  be  averaged  “physical” 
energy  (cf.  (2.30)) 


£- 


F-lW 


rt+1 


V 


n  +  1 


meas(f2^) 


(5.23) 


with  Ei  the  element  contribution  to  the  global  energy  E  and  meas(^)  the  area  of  element  L. 


Element  error  indicators  for  mesh  refinements  r\K  for  typical  time  step  are  chosen  in 
such  a  way  that 
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(5.24) 


1/2 

Kife  £  (X*  ) 

Accordingly,  %  are  computed  as: 

^  =  ^k.  2  +  ^2^k,  i 
+  ICJ^  )T]tc  2 

for  2- stage  scheme, 

_2  _  —2  ,  2  ^,2^2  2 

*7*  ”  *7*.  3  +  ^3^.2  +  ^*3^2^.  1 

+  2(q4'}  +  CjC^i0)^^ 

+  2C3C2£(11)rj/r2 

for  3-stage  scheme,  and 

_2  _2  ,  /.2J  ,  ,  r2n2r2J2 

Vk  ~  ^k. 4  +  ^4^.3  +  +  CjC-tC'flx  i 

+  2(C4£<1}  +  C4C3£<1)  +  C4C3C2£(11))^i4 

+  2(C42C3£<1)  + 

+  2C4C3C2£j1)  t)k  2 


(5.25) 


(5.26) 


(5.27) 


for  4-stage  scheme. 

The  solution  strategy  for  the  “one-stage  problem”  deserves  special  attention.  Since  the 
formal  operator  equivalent  to  the  bilinear  form  B\  has  zero  off-diagonal  terms,  velocity 
equations  and  displacement  equations  can  be  solved  independently  and  in  parallel  (the 
coupling  between  velocity  and  displacements  is  only  through  the  right-hand  side). 
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Finally,  we  emphasize  the  need  of  using  an  iterative  linear  equation  solver  in  the 
adaptivity  loop  (cf.  Step  12).  Notice  that  the  matrix  defining  the  left-hand  side  of  the 
resulting  algebraic  problem  is  syrr  netric  and  positive  definite  which  is  the  most  popular 
case  in  the  theory  of  iterative  solvers.  Notice  also  that  the  total  computational  effort 
involved  in  the  adaptivity  loop  is  governed  by  the  total  number  of  iterations  required  to 
solve  the  resulting  system  of  linear  equations  on  different  meshes  rather  that  by  the  total 
number  of  mesh  adaptations. 


6  Numerical  Example 


Stress  wave  propagation  in  an  elastic  panel  with  a  cut-off  (2D) 

As  an  example  we  chose  the  problem  of  wave  propagation  in  an  elastic  panel  with  a 
cut-off  (Fig.  13).  A  plane  progressive  transient  wave  is  normally  incident  on  a  traction-free 
rectangular  cut-off  and  gives  rise  to  complicated  pattern  of  stress  field  caused  by  the 
interference  of  incident,  reflected  and  diffracted  waves,  and  by  singularities  at  the  comers. 

The  problem  was  solved  for  the  following  data: 

•  Lame  coefficients,  p  ~  0.25  Pa,  A  =  0.5  Pa 

•  Mass  density,  pQ  =  1  kg/m 3 

•  Initial  conditions : 


jc-velocity,  vx(x,  y,  0)  ~f(x)[H(x  -  1)  -  H(x  - 1 )] 
y-velocity,  v^ft,  y,  0)  =0 

x-displacement,  ux(x,  y,  0)  =  F(x) 
y-displacement,  uy(x,  y,  0)  =  0 


(6.1) 


where  ffx)  =  - 1  +  32(|  -  x)2  -  256(|  -  .t)4  and  F(x)  =  anti-derivative  offfx). 

The  problem  was  solved  by  using  a  2-stagc  4th  order  Taylor-Galerkin  method  with 
stability  parameter  A  =  0.47125  and  constant  time  step  At  =  0.0078125  s.  No  local  p- 
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refinements  were  used  in  this  example,  as  the  p  capabilities  of  the  code  were  used  only  to 
specify  an  initial  order  of  approximation  (p  =  4  in  the  example). 

Figure  14  shows  the  initial  finite  element  mesh  and  Figs.  15  and  16  show  the  initial 
condition  functions  in  the  form  of  3D  plots  (;c- velocity  and  x-displacement,  respectively). 
Figure  17  presents  the  adapted  finite  element  mesh  at  time  t  =  0.25  s  and  Figs.  18  -  23 
present  the  corresponding  distributions  of  velocity  components  ( vx  and  vy )  and 
displacements  ( ux  and  uy).  Displacement  vector,  as  the  most  interesting,  is  shown  in  the 
form  of  both  the  3D  plots  and  contour  maps;  velocity  components  are  presented  in  the  form 
of  3D  plots  only.  As  can  be  seen,  a  doubling  of  the  amplitude  of  Jt-velocity  occures  at  the 
part  of  the  left  face  of  the  cut-off  which  is  in  perfect  agreement  with  the  plane  wave  solution 
(the  wavefiel  is  locally  plane  at  this  region  of  space-time).  Also,  sharp  gradients  of  vy  and 
uy  in  the  vicinity  of  point  A(;t  =  1.5  m,y  =  0.5  m )  can  be  noticed.  Similarly,  Fig.  24 
presents  the  finite  element  mesh  at  time  t  =  0.5  s  and  Figs.  25  -  30  present  corresponding 
distributions  of  vx,  Vy,  ux  and  uy. 

Figures  31-33  show  the  finite  element  mesh  at  time  t  -  0.75  s  and  Figs.  34  -  40 
show  the  corresponding  solution  (vx,  vy,  ux  and  uy).  Figure  32  presents  a  zoom  of  the 
finite  element  mesh  in  the  vicinity  of  comer  A.  A  total  number  of  19  mesh  refinements  was 
needed  to  capture  the  singularity  at  A.  Thus,  the  CFL  number,  v : 


(6.2) 


is  of  the  order  of  103.  This  would  place  a  severe  limitation  on  the  time-step  size  for  any 
conditionally  stable  scheme.  Here  hK  is  the  element  size,  pK  the  corresponding  order  of 
approximation,  and  XK,x,  XK^y  are  the  maximum  values  of  eigenvalues  for  Jacobian 
matrices  Rkl.  Similarly,  Fig.  33  presents  a  zoom  of  the  finite  element  mesh  in  the  vicinity 
of  point  B(;t  =  1.75  m,y  =  0.5  m). 

Figures  40  -  42  show  the  distribution  of  the  components  of  stress  tensor  Txx,  Tyy  and 
Try,  in  the  form  of  3D  plots.  A  very  strong  singularity  in  stresses  at  the  comer  points  A 
and  B  are  observed.  Finally,  Fig.  43  shows  the  history  of  the  error  bound  (5.22)  for  the 
whole  evolution  problem,  and  Fig.  44  shows  the  number  of  mesh  adaptations  as  a  function 
of  time. 
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APPENDIX:  Coefficients  for  TG  Schemes 

Listed  in  Table  A1  are  coefficients  fij  and  Vj  of  factorization  (3.32)  for  various  schemes  and 
various  values  of  stability  parameter  X. 


Table  A1 


J 

m 

X 

vi  _  _  . 

2 

4 

0.47125 

-0.429094784839126 

1.42909478483913 

-0.378098940635512 

0.548816059850774 

0.956435464587639 

-1.24942270597577 

2.24942270597577 

0. 

1.39760887500830 

3 

6 

0.695 

1.49387844844708 

1.49654973431122 

-1.99042818275830 

0.646649962127165 

0.455862817625521 

1.02905635860739 

1.41218087134444 

1.59760032454113 

2.26129033257844 

-2.85889065711958 

0. 

1.55344122448889 

2.12952443875392 

4 

0.22 

1.20431944918114 

0.487902005145560 

-0.692221454326699 

0.423802122065587 

0. 

0. 

4 

8 

0.92 

1.51162577127872 

-0.834510555244671 

-2.07038751330515 

2.39327229727110 

0.699029877978738 

-0.525519943835058 

1.17729970655410 

1.46702890097419 

1.86773692752800 

2.24345154354026 

-1.80258065404047 

-2.9580687867851 1 

3.51719789728532 

1.63319393127839 

0. 

2.41911998051307 

3.17835029737392 
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6 

0.29 

1.40623086480151 

-1.45772045255471 

0. 

1.05148958775320 

0.492263592133862 

0.562362275582733 

0. 

0.389406390318465 

5 

10 

1.15 

1.82704798577991 

-2.74777882789812 

-2.09945299083185 

2.49966934516921 

1.52051448778085 

0.581806316107477 

1.91908830024251 

1.27190423205347 

1.66880094503773 

0.736603622974089 

2.32321418507049 

2.22324138974647 

-3.93321086941255 

-2.98201123343815 

3.63234175933788 

2.05963895376635 

1.68098156751573 

3.93536053161087 

2.58923026751565 

3.54903685460515 

0. 

8 

0.38 

0.641021237097618 

-0.905360775899868 

-1.47520238954498 

1.70120453471140 

1.03833739363583 

0. 

0. 

0.625291329822029 

0.751557791280019 

0.412675810404174 
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4,  X  =  0.47048. 


Fig.  2  Eigenvalues  of  transient  operator  T,  s  =  3,  m  =  6,  X  =  0.695. 


f  transient  operator  T,  s 


Fig.  6  Eigenvalues  of  transient  operator  7\  j  =  4,  m  =  6,  X  =  0.29. 
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Fig.  13  Wave  propagation  in  an  elastic  panel.  Problem  definition. 
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Fig.  14  Wave  propagation  in  an  elastic  panel  problem.  Initial  FE  mesh. 
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Fig.  16  Wave  propagation  in  an  elastic  panel  problem.  Initial  condition, 
the  x  component  of  the  displacement  vector. 
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Fig.  18  Wave  propagation  in  an  elastic  panel  problem.  The  x  component 
of  the  velocity  vector  at  time  t  -  0.25  s. 
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Fig.  21  Wave  propagation  in  an  elastic  panel  problem.  The  y  component 
of  the  displacement  vector  at  time  /  =  0.25  s. 


454R-04 


10192 


PROJECT:  elastl2 


mrasni 

isss8:c:s:s2B8s»®iisss^»&si 
igigw»'~»8aSSB88«38«^^^l 
■|4  Ssssss  ibkKssssii 


gg!$£5? 


WM  fiHMDW  MHkJU  >■ . 

Sm  &££&  tet  wfarae  I 

fgggtjp  ! 


adap 


PROJECT:  clast  1 2  FIRST  COMPONENT 


r- 

CN  w-j 
On  m 
oo  un 
<N  ON 


3S 

^  o 


Fig.  25  Wave  propagation  in  an  elastic  panel  problem.  The  a:  component 
of  the  velocity  vector  at  time  t  =  0.5  s. 
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Fig.  26  Wave  propagation  in  an  elastic  panel  problem.  The  y  component 
of  the  velocity  vector  at  time  t  =  0.5  j. 
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Fig.  27  Wave  propagation  in  an  elastic  panel  problem.  The  x  component 
of  the  displacement  vector  at  time  t  =  0.5  s. 
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Fig.  33  Wave  propagation  in  an  elastic  panel  problem.  An  h  adaptive 
finite  element  mesh  at  time  t  =  0.75  s.  Zoom  of  Fig.  31. 
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Fig.  35  Wave  propagation  in  an  elastic  panel  problem.  The  y  component 
of  the  velocity  vector  at  time  t  =  0.75  s. 
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Fig.  36  Wave  propagation  in  an  elastic  panel  problem.  The  *  component 
of  the  displacement  vector  at  time  t  =  0.75  s. 
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Fig.  37  Wave  propagation  in  an  elastic  panel  problem.  The  y  component 
of  the  displacement  vector  at  time  /  =  0.75  s. 
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Fig.  38  Wave  propagation  in  an  elastic  panel  problem,  ^-displacement 
contours  at  time  t  -  0.75  s. 
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Fig.  39  Wave  propagation  in  an  elastic  panel  problem,  ^-displacement 
contours  at  time  t  ~  0.75 
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Fig.  42  Wave  propagation  in  an  elastic  panel  problem.  The  x  - 
component  of  the  stress  tensor  at  time  t  =  0.75  s. 
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Fig.  43  Wave  propagation  in  an  elastic  panel  problem.  The  bound  on  the 
global  energy  norm  of  the  error  as  a  function  of  time. 


Fig.  44  Wave  propagation  in  an  elastic  panel  problem.  The  number  of 
mesh  adaptations  as  a  function  of  time. 


